#include "uhm.hxx"

#define N 1

using namespace uhm;

typedef ID_Node_H_<N>          Id_Nd;
typedef ID_Element_<N>         Id_El;

#ifdef NOPIV
typedef LU_nopiv_Linal_        Uhm;
#endif

#ifdef PIV
typedef LU_piv_Linal_          Uhm;
#endif

#ifdef INCPIV
typedef LU_incpiv_Linal_       Uhm;
#endif

typedef Direct_Solver_<Id_Nd,Id_El,Uhm> Solver;

int main (int argc, char **argv) {

  initialize(argc, argv);

  // ** Environments
  Int_
    type = Mat_Base_::DOUBLE,
    uplo = Mat_Base_::FULL,
    sym  = Mat_Base_::UNSYM,
    rhs  = 1,
    blk  = 256;

  Int_
    kind = Mesh_Base_::NODE_DEFAULT;

  // ** Solver initialization
  Solver solver(type, uplo, sym, rhs, blk);

  // ** Import a mesh from a file_in
  option_begin();
  std::string fin;
  get_option_string("-in", "Reading a file", __FILE__,
                    "../../../../uhmfile/toy_1.uhm", fin);
  Double_ err;
  get_option_double("-err", "Error level to determine PASS/FAIL", __FILE__,
                    1.0e-3, err);

  option_end();

  time_begin();

  std::cout << argv[0] << "  Opening file ... " << fin << std::endl;
  std::ifstream in;
  in.open(fin.c_str());
  stream2mesh<Solver::_Mesh_>(in, solver.mesh());
  in.close();

  // ** Create random matrices in work
  UHM_ERROR((mesh2cache<Solver::_Mesh_,Solver::_Cache_>(solver.mesh(), solver.work())),
            ">> mesh2cache Error");
  
  solver.work().set_matrices(type, rhs, blk);
  solver.work().create_matrices();
  solver.work().random_matrices(true);

  // ** create tree and allocation
  time_in("direct:: Creating Solver (analysis)");
  UHM_ERROR(solver.create(),    ">> Error in Creating");
  time_out();

  time_in("direct:: Decomposing");
  UHM_ERROR(solver.decompose(), ">> Error in Decomposing");
  double t = time_out();

  Double_ flop, recorded;
  solver.get_flop(flop, recorded);

  std::cout << "\n\nFLOP = " << (flop/1.0e9)
            << ", Updated  = " << (recorded/1.0e9) << "  GFLOP\n";
  
  std::cout << "\n\nFLOPS= " << (flop/1.0e9)/t
            << ", Modified = " << ((flop-recorded)/1.0e9)/t << "  GFLOPS\n";

  Double_ storage, max_storage;
  storage_counter(storage, max_storage);
  std::cout << "\n\nStorage = " << (storage/1.0e6) 
            << ", " << (max_storage/1.0e6) << "  MBytes\n\n";


  time_in("direct:: Solving");
  UHM_ERROR(solver.solve(), ">> Error in Solving");
  time_out();

  time_in("direct:: Checking");
  Double_ res;
  UHM_ERROR(solver.check(res), ">> Error in Checking");
  time_out();

  std::cout << "\n\nResidual = " << res << "\n\n";
  Int_ r_val = (res < err ? 0: -1);
  std::cout << (res < err ? "[PASS]\n\n" : "[FAIL]\n\n") ;

  UHM_ERROR(solver.clear(), ">> Error in Clearing");
  UHM_ERROR(solver.flush(), ">> Error in Flushing");

  time_end();

  finalize();

  return r_val;
}


